% =========================================================================
%% 
% COMPARE_IRF
%
% Simple version to compare IRFs of the LMJ project. 
% Allows for: 
% 1) Single mode and a perturbation of a single parameter, given user
% defined values 
%
% 2) Multiple modes 
% 
% 3) Single mode with no perturbation 
%
% If loading a single mode, then must provide a paremeter name and a range
% of values to perturb it. If left empty it will compute a single set of
% IRFS. 
%
%% Documentation User Defined Settings 
% 
%% A. File and model settings 
% 
% *LOADPATH:* Subfolder within this path from where will pick up XLS file
% with the mode or modes. E.g. 'Gnp\L2YChi1Gtv\dtrsep2010\09 13 10'
% 
% *OUTPATH_NAME:* Name of the output path, which will be a subfolder of
% LOADPATH.  E.g. 'T2d Perturb'
%
% *FNAME:* name of the XLS file (extension not needed). E.g. 'Modes'. 
%
% *LEGENDS:*, name of the different modes to be loaded (applies only to
% multiple modes), leave empty if loading a single mode. E.g. {'Mode
% 1','Mode 2'}; 
%
% *SHEETNAME:*, (optional), sheet name. Leave empty if using a single sheet XLS
% file, specify otherwise. 
%
% *FUNCMOD:*, function handle of which model to use for the solution, e.g.
% @MODGNPL2Y or @MODGNPL2YCHI1.  
%
%% B. Settings for the IRFS 
% 
% *IRF_REPORT:* cell with the name of the state series that will be reported
% in the IRF. See STANAMES within the model code. 
% *NOTE:* in order to obtain an IRF for the level of a series
% written as a growth rate in the code type its name following by lev, in a single cell. 
% E.g. for series 'yobs' (observable output) choosing 'yobs lev' will automatically do a cumsum of the IRFS. 
% 
% *IRF_SHOCKS:* cell with the name of the shocks for which to report 
% the IRFS. If left empty, all shocks will be reported. E.g. ={'Neutral
% sta','MEI sta'}; 
%
% *IRF_REPORT_GNAMES* (optional) cell with alternative names wish to use for the
%  graphs for the states defined in *IRF_REPORT*. Both must be of the same length.
%  E.g. if IRF_REPORT={'yobs lev',u'}, could have
%  IRF_REPORT_GNAMES={'Output','Unemployment'}, and the latter will be used in
%  the graphs. 
%
% *NIRFPER:* scalar with the number of periods for which to plot the IRFS.
% E.g. 24. 
%
%% C. Settings for single parameter perturbations 
%
% Note: if single mode loaded and *PERTURB_NAME* empty, [], then only IRFS to the
% mode will be computed. 
%
%
% These inputs apply only to the case in which a single parameter is
% loaded and one of its elements will be perturbed. If these are defined
% and a multiple parameter spreadsheet is loaded, then these settings are
% irrelevant, but will cause no error message. 
%
% *PERTURB_NAME*, string or cell with the name of the parameter that will be perturbed 
% E.g. 'T2d'. 
%
% *PERTURB_GRID*, vector with parameter values for the parameter specified
% in PERTURB_NAME. *Note:* the value of the mode will be included in the grid
% by default. Values will be sorted in ascending order, including the mode.
% E.g. =[7,20]; 
% 
%% D. Additional settings for the plots 
% 
% The following give great control over how the multiple IRFS will be
% plotted. If not defined or left empty, default values are assigned. 
% *See COMPARE_IRF_PLOTS.m* for details. 
%
% *IN.TITFONT:*Size of title  fonts 
%
% *IN.LEGFONT:* Size of legend fonts
%
% *IN.FLAG_ADDPATH:*, if ==1 will create a cover sheet with the path where the IRFS
% stored. 
%
% *IN.ADDCOVER:*, Cell (as large as 4x1) with text to put on the cover
% sheet 
% 
% *IN.FLAG_LEGENDOFF:* ==1, no legend will be produced, otherwise the
%    cell in *legends* will be put at the bottom of the graph
%
% 
%%  Notes
% Created Sep 13 2010 
%
% Alejandro Justiniano 
%
% Modifications: 
%% Settings to allow for graphs for paper/presentations Sep 17 2010 
%
% 1) If single mode loaded and *PERTURB_NAME* empty, [], then only IRFS to the
%    mode will be computed. 
%
% 2) If *IN.FLAG_LEGENDOFF=1* No legend will be produced, otherwise the
%    string in *legends* will be put at the bottom of the graph 
%
% 3) *IRF_REPORT_GNAMES* (optional) cell with alternative names wish to use for the
%    graphs for the states defined in *IRF_REPORT*. Both must be of the same length. 
%    E.g. if IRF_REPORT={'yobs lev',u'}, could have
%    IRF_REPORT_NAMES={'Output','Unemployment'}, and the latter will be used in
%    the graphs. 

%
%% For Claudio
%
%% Code names These are the names of the main codes currently at use 
%
% *1) mod   Gnp  L2Y       Gtv* 
%
% *2) mod   Gnp  L2Y Chi1   Gtv* 
%
% L2Y :  parameter is directly on L2Y100, RBAR is implicit 
% Chi1:  Chi==1, first code, Chi > 1 
% Gtv :  The G solution matrix is time varying, since all variables are
% monthly and use the cummulator representation for the mixed frequency
% estimation. 
%
%% Main files
% 
% The current prior for botjh model is *dtrsep2010 Sh'. 
% 
% The name of the main output folders are as follows: 
% 
% For *modGnpL2YGtv* 
%
% *\Gnp\L2YGtv\dtrsep2010\09 10 10:* this is the long-run over the
% weekend. See tab_loop.xls and complete output there. 
% The refined run of that is *Gnp\L2YGtv\dtrsep2010\09 13 10*. See
% *tab_loop.xls*. The 3 modes from this run are stored in *modesSep13.xls*.
% 
% 
% For *modGnpL2YChi1Gtv:* 
%
% *Gnp\L2YChi1Gtv\dtrsep2010\09 10 10:* this is the long-run over the
% weekend. See tab_loop.xls The refined run of that is
% *Gnp\L2YChi1Gtv\dtrsep2010\09 14 10*, see tab_loop.xls. 
% The 2 modes from this run are stored in *modesSep14.xls*. 
% 
%
%
% For the names of states and shocks, refer to outputs *STANAMES* and
% *SHONAMES* withing each code. 
%
clear all; close all; cucd=cd; 
% =========================================================================
%% USER SETTINGS 
%
% Unless otherwise explained here, refer to descriptions above 
% =========================================================================

%% A.i) Files and names 
cucd = pwd;

loadpath='O:\PROJ_LIB\LJM\exog\US\MEAll\Cubic NormalRatio 10\Prior on AvWithin 83 Restart 11 02 11'; 
funcmod =@modUSMEAll; 
outpath_name ='DebugMatt'; 
%fname='modes'; % xls file that stores the modes
fname='TwoModesDeblanked'; % xls file that stores the modes
set_fname = 'settings';  
%autoset = 1; 
% legends={'Higher','Lower'}; 
legends={'Tight','Cal2'}; 
sheetname=[]; 

outfname = 'analysis'; 


%% A.ii) Model and solution settings 
% Call settings file
cd(loadpath);
eval(set_fname);
cd(cucd);

%% load in data; 
cd(loadpath); 
load workspace Y topcell 
cd(cucd); 
if ~exist('topcell','var'); topcell = []; end; 

%% B) Settings for the IRFS 
%irf_report
%if autoset == 0;

%% C) Settings for single parameter perturbations 
% redundant if multiple modes loaded 
perturb_name   =[]; 
perturb_grid   =[];  % Loaded value will be included 
%% D) Settings for plots
in.titfont=12; 
in.legfont=12; 
in.flag_legendoff=1; 
in.flag_addpath=1; 
% in.addcover={'Highest mode Chi =1 GTV';'Second mode very different wage rigidity'}; 
in.addcover={'tight';'cal2'}; 
% =======================================================================
% ==================== END USER INPUTS ===================================

%% I. Load xls file FNAME from sheet sheetname

disp(['Name of model:',func2str(funcmod)]); 
pause(0.5); 

cd(loadpath);
if exist('sheetname','var')==0 || isempty(sheetname)
    [parmodes,parnames]=xlsread(fname);
else
    [parmodes,parnames]=xlsread(fname,sheetname);
end
if size(parmodes,2)==1
    disp('Loaded single parameter value');
    pause(0.2);
    flag_single=1; 
else
    disp('Loaded multiple parameter values');
    pause(0.2);
    flag_single=0; 
end
cd(cucd);
printcell([parnames(:) num2cprec(parmodes)]);  
outpath=cr_dir(loadpath,[outpath_name,' ',strdate] );



%% II. Check model solves at first column of loaded values 
% Use this to get stanames and ssnames 

[GG, RR, CONS, eu, SDX, ZZ, outfield, ssvec, flag, ssnames, stanames,shonames]...
     = feval(funcmod,parmodes(:,1), solveopt, addsol);
 if ~isequal( eu,[1;1] ) 
     error('Cannot solve model at first column of parmode') 
 end 
 if length(shonames)~=size(SDX,1) 
     error('Wrong dimension of shonames') 
 end 
 
 %% Check IRF variables and shock names exist 
 %  Also check IRF_REPORT_GNAMES 
 disp('Check states to be reported exist');
 state_pos=extractlevels(state_sel,stanames);
 disp('Check shocks to be reported exist');
 if exist('irf_shocks','var')==0 || isempty(irf_shocks); 
     shocks_pos=(1:size(SDX,1)); 
 else 
     shocks_pos=cellposition(irf_shocks,shonames); 
 end 
 
 if exist('state_sel_gnames','var') && ~isempty(state_sel_gnames);
     disp('Will use alternative name for IRFS');
     disp('Original Names      Graph Names');
     if length( state_sel(:)  ) ~= length( state_sel_gnames(:) )
         dispaj('Length IRF_REPORT:',length( state_sel(:)  ) );
         dispaj('Length IRF_REPORT_GNAMES:',length( state_sel_gnames(:)  ) );

         error('Lengths of IRF_REPORT and IRF_REPORT_GNAMES do not match')
     end
     printcell([state_sel(:) state_sel_gnames(:)]);
     pause(1);
     
 end
 
 
%% II. If using a single parameter perturbation, create grid in ascending
%order 

if flag_single
    if ~isempty( perturb_name )
        disp(['Perturbing :',perturb_name]);
        perturb_pos=cellposition(perturb_name,parnames);
        perturb_grid=sort([perturb_grid(:);parmodes(perturb_pos)],'ascend');
        parmodes=repmat(parmodes,[1 (length(perturb_grid))]);
        parmodes(perturb_pos,:)=perturb_grid;
        legends=fnumcell([perturb_name,' '],perturb_grid,3);
        disp('Perturbation values');
        printcell( num2cprec(perturb_grid(:)) );
    else
        display('Perturbname is empty, perturbing a single mode')
        legends={'Mode'}; 
    end
end   

%Create PDF of Parameters
cd(outpath)
if isfield(outfield,'parTex')
    parTex = tex2string(outfield.parTex,parnames);
else
    parTex = parnames;
end
MakeLaTeX_table(parmodes,parTex,legends,'','Parameter Values','paramTable',0)
cd(cucd)

%% III. Solve model for each value of PARMODES, also store SS vector 
nstates=length(ssvec); 
ncols    =size(parmodes,2); 
nirf     =length( state_sel ); 
nx       =size(SDX,1); 

ss_mat=zeros(nstates,ncols); 
irf_mat=zeros(nirf,nirfper,nx,ncols); 
ok_vec=zeros(ncols,1); 

for ii=1:ncols;
    [junk,junk,irfstates,ssvec]=irf_quick(funcmod,parmodes(:,ii),solveopt,addsol,nirfper,state_sel,0);
    if ~isempty(irfstates)
        ok_vec(ii)=1;
        irf_mat(:,:,:,ii)=irfstates;
        ss_mat(:,   ii)=ssvec;
    end
end

% Delete those for which could not find a solution 
pos_ok=find(ok_vec==1); 
irf_mat=irf_mat(:,:,:,pos_ok); 
ss_mat =ss_mat(:,pos_ok); 
legends=legends(pos_ok); 


%% IV. Create IRFS, dps file of IRFS 
ghvec=compare_irf_plots(irf_mat,shocks_pos,state_sel_gnames,shonames,legends,outpath,in); 

%% V. Create XLS with SS values
cd(outpath); 
 sscell=crtcell( ss_mat, ssnames , legends ,{'SS Values'} ); 
% xlswrite('SS tab',sscell); 
if isfield(outfield,'ssTex')
    
    ssTex = tex2string(outfield.ssTex,ssnames);
    
else
    ssTex = ssnames;
end

MakeLaTeX_table(ss_mat,ssTex,legends,'','Steady States','ssTable',0)
cd(cucd); 

disp('SS for different modes'); 
printcell(sscell); 

disp('_______________________'); 
disp('COMPARE_IRF complete'); 
display(['Output saved in :',outpath]); 

close all; 

%% VI. Starting with other tasks
% that were to be found in estima_out. 
% The goal is to set up and run lmj_modeanalysis_compute


%% Spectrum

    
for ii = 1:ncols
    parmode = parmodes(:, ii);
    if flag_mixfreq == 0;
        [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames, stanames ,shonames]...
            =feval(funcmod,parmode,solveopt,addsol);
    else
        % for mixed frequency
        [GG, RR, CONS, eu, SDX, ZZ_zero, out, ssvec, flag, ssnames, stanames ,shonames]...
            =feval(funcmod,parmode,solveopt,addsol);
        ZZ = [ZZ_zero;  out.AA];
    end
    
    if ii == 1;
        znames=cell(size(ZZ,1),1);
        for zz=1:length(znames)
            temp      =find(ZZ(zz,:)==1 ) ;
            if length(temp) > 1
                error('Need 1 entry only for each row of ZZ')
            end
            znames(zz)=stanames( temp );
        end
    end
    
    figname = [legends{ii}, '_specdecom'];
    if flag_mixfreq == 0
        [out,spec_outcell,hand_vec]=dsge_spectrum(GG,RR,ZZ,SDX,specrep,stanames,shonames,znames,...
            ana_in.lowper,ana_in.highper,ana_in.npoints,outpath, figname);
    elseif flag_mixfreq == 1;
        if flag_gtv == 1;
            ZZ = [ZZ_zero;  out.AA];
        else
            ZZ = ZZ_zero;
        end
        [outq,outm,spec_outcell]=dsge_spectrum_mixfreq(GG,RR,ZZ,SDX,specrep,stanames,shonames,...
            znames,ana_in.lowper,ana_in.highper,ana_in.npoints,3,outpath,figname);
    end
    
    tableOpt.orientation = 'landscape';
    MakeLateX_table(out.vfbdecst,specrep,shonames,'','Variance Decomposition',...
                    ['varDec',legends{ii}],0,outpath,tableOpt);

    spec_outcell = merge_cells(topcell, spec_outcell, 2);
    close all
    cd(outpath); 
    try
        xlswrite(outfname,spec_outcell,[legends{ii}, '_spectrum']);
    catch
        disp('Spectrum Tab not saved')
    end
    cd(cucd); 
    
    
    %% Moments 
    
    % Theoretical 
    [vdcom,vcvth,sdevf,ACzero,ACone,vcshocks,auone,vdech,vdech_st,vdcom_st]=vardecom_mode(GG,RR,SDX,ZZ,50);
    var_outcell = crtcell( vdcom, znames , shonames ,['Stationary Variance_', legends{ii} ]);
    var_outcell = merge_cells(topcell, var_outcell, 2);

    cd(outpath); 
    try
        xlswrite(outfname, var_outcell,['stationary_', legends{ii}]);
    catch
        disp('Stationary Variance Tab not saved')
    end
    cd(cucd); 
end
delete  'variance decomposition*'

% =========================================================================
% Begin Simulated Moments 
%%

if ~exist('ana_in', 'var'); % ANA_IN does not exist as a variable
    ana_in = []; 
end
[e, ana_in] = ch_field(ana_in, 'flag_addz',  zeros(1,size(Y, 2)));
ana_in.flag_mixfreq = flag_mixfreq; 
ana_in.flag_gtv = flag_gtv; 
ana_in.stanames = stanames;
ana_in.shonames = shonames;
ana_in.parnames = parnames;
ana_in.state_sel = acstate_sel; 
%ana_in.modelnames = {'higher mode', 'lower mode'};
ana_in.modelnames = legends; 
ana_in.ssnames = ssnames;
ana_in.savepath = outpath;
ana_in.flag_acov = 0; 
ana_in.flag_acdecomp = 0; 
display('Cannot handle mixed frequency yet...'); 
ana_in.case_tag=[]; 
%ana_in.case_tag = outfield.case_tag; 
ana_in.nper = 4;   
ana_in.flag_medbonly = 0; 

if flag_mixfreq == 1
    [out_tag,cellrep]=tag_admissable(outfield.case_tag);
end
ana_out = lmj_modeanalysis_compute(funcmod, parmodes, Y, ana_in, solveopt, addsol);
pertmat = ana_out.pertmat;


% =========================================================================
% =========================================================================
%% Begin Output Section


if flag_mixfreq == 0; 
    stdsimz_mat = ana_out.stdsimz_mat; 
elseif flag_mixfreq == 1; 
    stdsimz_tag_mat = ana_out.stdsimz_tag_mat;
    stdsimz_mat = ana_out.stdsimz_mat;
end 

if ~ch_field(ana_out,'stdsimz_tag_mat')
    display('STDSIMZ_TAG_MAT does not exist.  Changing to []');
    stdsimz_tag_mat=[];
end

if  ~isempty(stdsimz_mat); % Base Frequency
    for ii = 1:ana_out.num_valid;
        stdsimz_ii = squeeze(stdsimz_mat(:, :, ii));
        std_outcell = playmod_outdraws(stdsimz_ii', znames);
        std_outcell = std_outcell(3:end, :); 
        std_data = ['Data'; num2cprec(std(Y))'];
        std_outcell = [std_outcell(:, 1), std_data, std_outcell(:, 2:end)]; 
        std_outcell = merge_cells({'Base Frequency (Monthly)'}, std_outcell); 
        if ~isempty(stdsimz_tag_mat) && flag_mixfreq == 1; 
            stdsimz_ii = squeeze(stdsimz_tag_mat(:, :, ii));
            std_outcell_tag = playmod_outdraws(stdsimz_ii', znames);
            std_outcell_tag = std_outcell_tag(3:end, :);
            std_data = ['Data'; num2cprec(std(tag_datasim(Y, out_tag)))'];
            std_outcell_tag = [std_outcell_tag(:, 1), std_data, std_outcell_tag(:, 2:end)];
            std_outcell_tag = merge_cells({'Aggregated Frequency (Quarterly)'}, std_outcell_tag);
            std_outcell = merge_cells(std_outcell, std_outcell_tag, 1);

        end
        
        if ii == 1;
            std_outcell_full = std_outcell;
        else
            std_outcell_full = merge_cells(std_outcell_full, std_outcell, 2 );
        end
    end
    std_outcell_full = merge_cells(topcell, std_outcell_full, 1);
    
    cd(outpath)
    xlswrite(outfname,std_outcell_full,'std_sim');
    cd(cucd)
    
    % 1) plot volatilities
    in_stdhist.modenames = ana_out.modelnames;
    in_stdhist.varnames  = znames;
    in_stdhist.savepath  = ana_in.savepath;
    in_stdhist.savename  = 'Z STD';
    in_stdhist.flag_ksd  = 1;
    plot_stdhist(stdsimz_mat, Y, in_stdhist);
    
end

% 3) plot Variance Decompositions.
% 3.1) of the observables
if ~isempty(ana_out.vdechz_mat)
    in_vdec.modenames = ana_out.modelnames;
    in_vdec.varnames =  znames;
    in_vdec.shonames = shonames;
    in_vdec.savepath = outpath;
    in_vdec.savename = 'Z VDEC';
    plot_vdecs(ana_out.vdechz_mat, in_vdec);
end

% 3.2) of the variables
if ~isempty(ana_out.vdechs_mat)
    in_vdec.modenames = ana_out.modelnames;
    in_vdec.varnames =  ana_in.state_sel;
    in_vdec.shonames = shonames;
    in_vdec.savepath = outpath;
    in_vdec.savename = 'Sel State VDEC';
    plot_vdecs(ana_out.vdechs_mat, in_vdec);
end

%%

% 4) AUTOCORRELATION
% 4.1) of the obervables;
in_ac.modenames = ['data';  ana_out.modelnames(:)];
in_ac.varnames = znames;
in_ac.savepath = outpath;
in_ac.savename = 'Obs AC Median';

if flag_mixfreq == 1
    [out_tag,cellrep]=tag_admissable(outfield.case_tag);
    ac_data=accov(tag_datasim(Y, out_tag),(0:ana_in.nper));
elseif flag_mixfreq == 0
    ac_data=accov(Y,(0:ana_in.nper));
end

if flag_mixfreq == 0; 
    acsimz_mat = ana_out.acsimz_mat;
elseif flag_mixfreq == 1; 
    acsimz_mat = ana_out.acsimz_tag_mat;
end

if ~isempty(acsimz_mat);
    clear med_acsimz_mat;
    for ii = 1:ana_out.num_valid + 1;
        if ii == 1;
            med_acsimz_mat(:,:,:,ii)=ac_data;
        else
            med_acsimz_mat(:,:,:,ii)=median(acsimz_mat(:, :, :, :, ii-1),4);
        end
    end
    plotmultiple_acorr(med_acsimz_mat, in_ac);
end

% 4.2) of the selected states
in_ac.modenames = ana_out.modelnames;
in_ac.varnames  = ana_in.state_sel;
in_ac.savepath = outpath;
in_ac.savename = 'Sel States AC Median';
in_ac.flagblock = 1; 
if ~isempty(ana_out.acsims_mat)
    clear med_acsims_mat
    for ii = 1:ana_out.num_valid;
        med_acsims_mat(:,:,:,ii)=median(ana_out.acsims_mat(:, :, :, :, ii),4);
    end
    plotmultiple_acorr(med_acsims_mat, in_ac);
end

% 4.3) full cross-correlogram vs. data
in_acfull.modenames =  [ana_out.modelnames(:); 'data'];
in_acfull.varnames  = znames;
in_acfull.savepath  = outpath;
in_acfull.savename = 'Obs AC full';
%in_acfull.title = {'Obs AC full'};
if ~isempty(acsimz_mat);
    %plotmultiple_acall(acsimz_mat, ac_data, in_acfull);
    plot_all_autocorr(acsimz_mat, ac_data, in_acfull)
end

%Print Everything to One PDF
cd(outpath)
nModels = length(legends);
fileCell = cell(4+nModels,1);

fileCell{1} = 'paramTable';
fileCell{2} = 'ssTable';
fileCell{3} = 'Z STDDensity';
fileCell{4} = 'Obs AC Full';

for k = 1:nModels
    fileCell{4+k} = ['varDec',legends{k}];
end

mergePDF(fileCell,'allFigs',1);
cd(cucd)